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Abstract 

Systems with time delay play an important role in modeling of many physical and biological 
processes. In this paper we describe generic properties of systems with time delay, which are 
related to the appearance and stability of periodic solutions. In particular, we show that delay 
systems generically have families of periodic solutions, which are reappearing for infinitely many 
delay times. As delay increases, the solution families overlap leading to increasing coexistence of 
multiple stable as well as unstable solutions. We also consider stability issue of periodic solutions 
with large delay by explaining asymptotic properties of the spectrum of characteristic multipliers. 
We show that the spectrum of multipliers can be splitted into two parts: pseudo-continuous and 
strongly unstable. The pseudo-continuous part of the spectrum mediates destabilization of periodic 
solutions. 
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I. INTRODUCTION 



The dynamical behavior of various physical and biological systems under the influence 
of delayed feedback or coupling can be modeled by including terms with delayed arguments 
in the equations of motion. When the delay becomes longer than the other characteristic 



timescales of the system, a complicated and high-dimensional dynamics can appear 
y, The analysis and control of such dynamical regimes is important for many applications 
including lasers with optical feedback or coupling Q, I^, neural activity control 0, [lo| . 
and many others [uj. For instance, the following complicated regimes have been observed 
in lasers with delayed feedback: low- frequency fluctuations [l^, regular pulse packages jl^, 
coherence collapse, just to mention a few. 

One of the fundamental question in the analysis of systems with delay concerns properties 
of periodic solutions. Periodic solutions of any dynamical system, including also systems 
with delay, are important parts of the dynamics. When such solutions are stable, they 
can be directly observed experimentally or numerically. In the case, when such solutions are 
unstable, they play an important role, e.g. by determining of a set of admissible initial values 
to be attracted to some stable steady state (basin boundary), or by forming fundamental 
blocks of a chaotic attractor 1^, [l^. Finally, unstable as well as stable periodic solutions 
can play an important role in mediating any kind of soft or hard transitions as some control 
parameters are varied. 

This paper is devoted to the study of generic properties of periodic solutions in systems 
with a constant time delay 

x'(t) = /(x(t),x(t-r)), (1) 

where x G i?" and r > is the time delay. Since we investigate the influence of the 
delay, we assume r to be our control parameter. Delay has been used as a parameter in 



various applications: chaotic systems with feedback 
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networks or arrays of oscillators with delaved coupling 
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17 1 , network motifs 
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, l26|, |27|, laser systems with feedback [28|| and couplin 



10, 



33, 
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HQ, 



large 



23l |. mechanical systems 
29 1 , coupled neurons [3( 
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39|, 



chemical oscillators [32i|, delayed feedback control [8|, 
believe that our results are applicable to the all above mentioned well as to many 

others, which include time delay as a controllable parameter. 

The plan of the paper is as follows. Section [Til starts by showing that periodic solutions 
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in system ([TD are forming branches with respect to the control parameter r. These branches 
are reappearing infinitely many times for different delay values. As the delay increases, 
the solution branches overlap leading to increasing coexistence of multiple stable as well 
as unstable periodic solutions. The number of coexisting solutions is shown to be linearly 
proportional to the delay time. Further in Section IIIIl we consider stability properties of 
periodic solutions for systems with large delay by explaining asymptotic properties of the 
spectrum of their characteristic multipliers. The spectrum of multipliers can be splitted into 



e. Such situation is similarl to the case of 



42| . The pseudo-continuous spectrum 



two parts: pseudo-continuous and strongly unstab 
steady states in delay systems with large delay js, 
mediates bifurcations of periodic solutions for systems with large delay. Sections HV] and |V] 
discuss some implications of the existence of pseudo-continuous spectrum and possibility for 
its numerical computation. 

The obtained results provide a better understanding of mechanisms behind the growing 
multistability and dynamic complexity in systems with delay. In particular, we show that 
coexistence of multiple stable (as well as unstable) periodic solutions is a natural feature 
of delay systems. The growing "effective dimension" of dynamics with the growing delay is 
supported by the fact that the dimensionality of unstable manifolds of periodic solutions 
also grows linearly with delay. Our results will be illustrated using Duffing oscillator with 
delay 

x"{t) + dx'it) + ax{t) + x^it) + h {x{t) - x{t - r)) = 0, (2) 
where x{t) is a real variable, and Stuart-Landau oscillator with delayed feedback 

z'{t) = (a + ip)z{t) - z{t)\z{t)\^ + Kz{t - r), (3) 

where z{t) is a complex variable. 



II. REAPPEARANCE OF PERIODIC SOLUTIONS AND BRANCHES 

In this section we show that periodic solutions of system ^ reappear infinitely many 
times for different and well specified values of delay r. We will also show that such solu- 
tions typically form branches, which can be mapped one onto another by some similarity 
transformation. 
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Figure 1: Reappearance of a periodic solution XQ{t) for the delays r„ = tq + nTg, n = 1,2,3, .. .. 
To is the period of xo(t). 

A. Reappearance of periodic solutions 

Let us consider system ([T]), which possesses a periodic solution Xo(t) for a time delay 
T = tq. Let To be the period of this solution. Then it is easy to check that the same periodic 
solution exists in system ^ with time delay ri = ro + Tq. Indeed, substituting xo{t) into 
([T]) we obtain 

x'oit) = f{xo{t), Xo{t - Ti)) = f{xo{t),Xo{t - To - Tq)) 

= /(xo(t),a;o(t - To)), 

where the periodicity of a;o(t) is taken into account. Similarly, the solution Xo(t) reappears 
for all values t„ = tq + nTo, n = 1, 2, 3, ... of the delay; see Fig. [H 

B. Reappearance of branches of periodic solutions 

A periodic solution Xoit) can be generically continued to a branch of periodic solutions 
Xo{t; t) with respect to the parameter r, at least in some neighborhood of To. Denote T(r) 
to be the period of these solutions along the branch. Since each individual periodic solution 
reappears for every delay time T + nT{T), the whole branch will appear infinitely many times 
as well; see Fig. [U It is naturally to introduce the notion of primary branch, which satisfies 
r < T(r). For convenience, let us parametrize the primary branch using the parameter /, 
which coincides with the delay on this branch Xo(t; /) := Xo(t; r). A solution Xo(t; I), which 
corresponds to some parameter value Z, will appear again on the n-th branch at time delay 
r(n, /) = / + nT{l); see Fig. [H Thus, we obtain the representation of the n-th branch 



Xn{t; T{n, I)) = Xn(t; I + nT{l)) = Xo(t; I). 



(4) 
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Figure 2: Branches of periodic solutions for Duffing oscillator. Left panel shows the amplitude 
of the solutions versus time delay, right panel shows dependence of the period along the branch. 
Parameter values: (a,b) a = 0.5, b = 0.6, d = 0.06; (c,d) a = 1.38, b = 0.4, d = 0.3. The first three 
branches are plotted in black and the rest in gray. 



The corresponding mapping, which maps delay times, is given as follows 



/^r(n,/) = l + nT{l). 



(5) 
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le above described branches can be numerously found in the research 
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Usually, these branches can be found numerically. 



A useful tools for finding such branches is the continuation software DDE-Biftool [4j|. In 
fact, only the primary branch should be computed while the others can be obtained using 
the transformation (Ilj). 

For the model of Duffing oscillator (JS]), we have found two types of branches, which are 
presented in Fig. [2l Other examples are given in Fig. [H Note that due to the nontrivial 
dependence of the period T(r) along the branch (see the right panel of the figure), the 
mapping ((51) is not just a parallel shift. It has some more complicated properties, which will 
be studied in the following subsection. 



C. Properties of the branches 

As we have seen, the primary branch of periodic solutions Xo{t;l) is characterized by a 
period function along the branch T(/). Clearly, since all other branches consist of the same 
solutions, they have the same period dependence T{t). Via the mapping ([5]), the function 
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T{1) determines uniquely how branches reappear for larger delay times. Let us discuss main 
properties of the map ([5]) and corresponding implications. 



Stretching, squeezing. 



Under the transformation 
squeezed. In particular if 



some parts of the branch will be stretched and some parts 
dT{n, I) 



dl 



\l + nT'{l)\ > 1 



(6) 



then the corresponding part will be locally stretched. Here T'(/) = dT{l)/dl is the derivative 
of the period function. If the inverse inequality 

dT{n, I) 



dl 



\l + nT'(l)\ < 1 



(7) 



is satisfied, the corresponding branch parts will be locally squeezed. With the increasing of 
n (which is equivalent to increasing r) almost all parts of the branches will be stretched, 
since ^ will be satisfied. Hence, the branches become eventually wider and occupy larger 
r intervals. This leads to the growing overlapping of branches and growing co-existence of 
periodic solutions with increasing r. This effect is clearly visible in Figs. [2](a) and [H 



Multistability. 

Let us describe the above mentioned phenomenon on a more quantitative level. We 
will show that the number of coexistent periodic solutions grows linearly with r and give 
estimations for the corresponding coefficient. Let us consider the two possible cases: 

CASE 1. The primary branch is confined to some interval of r as in the case shown in 
Fig. [2](a-d) . This means that the primary branch ranges from /min till /max (/min can be zero). 

CASE 2. The primary branch is bounded only from below, i.e. it ranges from /min till 
+oo; see an example in Fig. [4](a). 

Consider the first case. Let us denote Tmax and Tmin to be the maximum and the minimum 
of the period function T(/) on the interval /mm ^ / ^ /max- If ^max = oo then the problem 
can be reduced to the Case 2, since the next branch Xi(t; / + T(/)) = Xo{t; /) will be stretched 
up to r = oo by the mapping ([5]). Hence, Tmax and Tmin can be considered to be bounded. In 
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Figure 3: Illustration to the derivation procedure of the formula for the number of overlapping 
branches at time delay r. Shown is the delayed Duffing oscillator with parameter values a = 1.38, 
h = 0.4, d = 0.3, see also Fig. [2](c). The branches k and m are shown in black. More details are 
given in text. 

this case, all other branches are also bounded and exist for delay values rin, I) = I + nT{l), 
Imin < I < ^min- In particular, the n-th branch ranges from 

'rmin(n,/)= min r(n,/)= mm [l + nT{l)] 

^min _:^^^max i'min_i''_:i''m.ax 



till 



Tmaxin,!) = max T{n,l) = max [l + nT{l)]. 

^miii^^^^max ^ rain. — ^^^max 



For large enough n, the minimal and maximal bounds of the branches can be well approxi- 
mated as follows 



Tmin{n,l) = nmin 
T-max(n,/) = nmax 



n 



-+m 

n 



nTrr 



(8) 
(9) 



up to the terms of order one. Let us now find how many branches are overlapping for some 
sufficiently large delay value r. It is clear that these branches should satisfy the condition 

I) <T < Tm!,x{n,l). 

Let m be the least number of the branch, which exists at time delay r and k be the 
largest number, i.e. 

Tmax("i, I) and T^inik, I) ^ T 

up to terms of order one, see Fig. [3l Taking into account ([8]) and ([9]), we obtain 
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where the approximation sign means that the equality is satisfied up to the order one terms 
{k and m are large). All branches with numbers m < n < k exist at time delay r, see Fig. [3l 
Hence, the number of overlapping branches can be estimated as follows 

^min ^max — Tmin 

N^k — m = k — fc ™™ = k — = KiT, (10) 

max max 

where the coefficient for this growth is given by 

J- max-' min min max 

Expression (fTOll gives also a lower estimation for the number of coexisting periodic solutions 
of system 1^. Indeed, if the branches are folded like in Fig. [2](a) then one branch may lead 
to more than one periodic solution for some delay values. 

In a similar way, one can show that the number of coexisting branches in CASE 2 can be 
estimated as 

N^H2T = -^r. (12) 

min 

Finally note that there may exist few primary branches, which cannot be mapped one onto 
another by the transformation (111). In this case, each branch will reappear with the increasing 
delay. The growth rate n fin this case is given as a superposition of the corresponding rates. 

Turning points on branches. 

The branches may have turning points, which correspond to a fold bifurcation for the 
family of periodic solutions. The condition for the branch n to have a turning point can be 
written as follows 



dl 



l + nT'(/) = 0. (13) 



T\l) = -1. (14) 

With the incresing of branch number n, the fold point tends to some asymptotic value, 
which is independent on n and given by the condition T'(/) = 0. 



or, taking into account 



Equation (fT3l) can be rewritten as 
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D. Example: Stuart-Landau oscillator with delay 



In this paragraph we consider the following system 

z'{t) = (a + ip)z{t) - z{t)\z{t)\^ + z(t - r) (15) 

with the instantaneous part as the normal form for an oscillator close to supercritical 
Andronov-Hopf bifurcation. Such system is called sometimes Stuart-Landau oscillator 

n n 

[3|, 145||. The additional term z{t — r) accounts for a delayed feedback. Here z{t) is a 
complex variable, i.e. the system has essentially two components, which can be chosen as 
real and imaginary parts of z{t). 

Due to symmetry properties of this system, some of its periodic solutions can be found 
analytically in the form of rotating waves Ae'^'^K Substituting this rotating wave into (fTSl ). 
we obtain the equation 

iuj = {a + -A'^ + e-*"^", 
which leads to the following expressions for the amplitude A and frequency uj 



A = y/a + cos LOT. (16) 

u = P — sin UT, (17) 

For the purposes of this paper, let us rewrite (fTBll and (fTTll in terms of the amplitude A and 
the period T = iT^joj: 



A= ^a + cos(2'K^j. (18) 
P - sin (271^) 

The expressions (fTSll and (fTOll are invariant under the change r — > r + riT, which reflects 
the fact that solutions on different branches are identical. 

Although the period T(r) along the branch is given in the implicit way by ( flQl l. one can 
obtain an explicit parametric representation of the branches with respect to T and r. For 
this, we introduce an additional parameter ip = Svrr/T. With the help of this parameter, 
solutions of (fT9l) can be written as follows 

TW = IT^^ (20) 
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Figure 4: Branches of periodic solutions for the Stuart-Landau oscillator with delay, which are given 
analytically in (|20D - (|2TD: fa) = h (h) = 2. 



ip + 27rn 



(21) 



/3 — sin '?/' 

Now the branches can be easily plotted by varying parameter ip. For /3 = 1 the branches 
are unbounded, see Fig. [4](a), and for /3 = 2 they are bounded, see Fig. [4](b). Moreover, for 
/? = 2 neighboring branches are connected. 

As it is expected, the coexistence of multiple periodic solutions grows with the increasing 
of delay. The number of coexisting branches can be estimated using (fTOll and (fT2l) as 



N 



1 13+1 

— T = — : r 



27r 



in the case < /? < 1 and 



1 



T ■ T 

mm max 



TT 



for /5 > 1. Taking into account folding of the branches, the number of periodic solutions 
grows twice as fast with the increasing of the delay, i.e. with the rates (/3 + l)r/7r and 2r/7r 
respectively. 
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III. STABILITY OF THE REAPPEARING PERIODIC SOLUTIONS, LONG DE- 
LAY ISSUES 



A. Some elements of the stability theory for periodic solutions 

Let us introduce necessary notations and shortly remind basic elements of the stability 
theory for periodic solutions of ([T|). 

The linearization of ([T]) around some periodic solution p(t) with a period T has the 
following form 

eit) = Ait)at) + B{mt-T), (22) 

where A(t) = Dif{p{t),p{t—T)) and -B(t) = D2f{p(t),p{t—T)) areT- periodic matrices nxn. 
Here Di and D2 denote partial derivatives with respect to the first and the second argument 
respectively. Any solution of (l22|l with some initial function ip{t) can be represented as 
x{t; s, (f) = \E'(t; s)(p, where '^{t; s) is the evolution operator j46l|. The monodromy operator 
is introduced as the evolution operator evaluated at the period 

U = ^(T;0). 

of the periodic solution p(t) is determined by a countable set of characteristic 



Stability 

multipliers [46|, l47|| /x-,, j = 1, 2, . . ., which are defined by the spectrum off/. The correspond- 
ing characteristic exponents are given as = |;Ln/ij. A periodic solution is asymptotically 
stable if all its multipliers have modulus less than one, or, equivalently, real parts of all char- 
acteristic exponents are negative. A bifurcation occurs when a multiplier crosses unitary 
circle parameter change. 

Practically, characteristic multipliers and stability of a periodic solutions can be computed 
using DDE-Biftool software. 



B. Stability of periodic solutions versus delay 

Considering r as the control parameter, stability properties of periodic solutions change 
as r is varied. In general, it is a challenging problem to find their stability, especially for 
larger r. Figure [5] shows largest characteristic multipliers (with the largest modulus) for 
a periodic solution of delayed Duffing oscillator ^ for different delay times. Even though 
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Figure 5: Largest characteristic multipliers of a periodic solution of delayed Duffing oscillator for 
different values of delay. Different figures correspond to the same solution on different branches: 
(a) branch n = 2 (r = 2.2); (b) n = 20 (r = 23.7); (c) n = 80 (r = 91.6), (d) n = 140 (r = 159.5). 
Other parameters are a = 0.5, b = 0.6, d = 0.06. 

the solution is the same, its stability properties change as one moves from one branch to 
another. One can clearly observe that more and more multipliers come close to the unitary 
circle making the problem more and more degenerate and numerically stiff for larger delay. 

In the following we propose an analytical technique, which overcomes the appeared dif- 
ficulty and shows how to approximate the characteristic multipliers for larger delay values. 



In particular we show that the characteristic multipliers have simi 

n 

the eigenvalues of stationary states for systems with large delay p. 



ar pro perties to that of 



42|. 



The linearization of equation ([T]) around the solution x„(t; r) on n-th branch (this solution 
coincides with xo(t; I)) at time-delay t = I + nT{l) has the form 



where 



at) = A{t-i)m+B{t-i)at-T), 

Ait- 1) = DJ{xo{t;l),Xo{t-l;l)), 
B{t-l) = D2fixoit;l),Xoit-l;l)) 



(23) 



(24) 



are T-periodic matrices, which depend only on function / and a shape of the solution 
Xo{t] I). It is important that A and B do not depend on the branch number n and the time 



12 



delay r, at which the system is considered. This allows us to study stability properties of 
periodic solutions asymptotically with t oo. Namely, by increasing r (or, equivalently, 
branch number n) we are actually jumping from one branch to another by keeping the 
relative position / within the branch fixed, see Fig. [H As a result, we consider the same 
periodic solution Xo(t;l), which exists for different infinitely increasing delay values and 
we study stability properties of this solution as r increases. In the following, we consider 
r to be continuous parameter and then apply the obtained results to the countable set 
r(n, /) = / + nT{l), n = 0, 1, 2, ... of delay values. 



C. Pseudo-continuous spectrum 

Here we show that periodic solutions possess a family of characteristic multipliers, which 
have the following asymptotic representation 

1^1.1 = l + -l{uj) + o(\) (25) 



with the increasing of delay. Here u is a parameter along the family. Using the analogy to 



42| . we will call such spectrum 



the spectrum of eigenvalues for stationary solutions [3|, w}. 
pseudo- continuous. Its main features are: 

1. Pseudo-continuous spectrum tends to the critical value = 1 as r — >^ cx); 

2. Its stability is determined by the sign of the function 7(^0'); 

3. For any finite r, the parameter u admits a discrete countable number of values. As 
T 00, the spectrum tends to a continuous in the sense that the discrete parameter 
uj covers densely the whole interval [0,27r]. 

In the case, if a periodic solution has only pseudo-continuous spectrum, its stability for large 
enough r will be uniquely defined by the function 7(0;). 

Before we give a rigorous proof for the existence of the pseudo-continuous spectrum, 
let us illustrate it using our numerical example of Duffing oscillator ([2]). Figure [6] shows 
approximations for the function 7(0;) by plotting 7^ = r(|/ij| — 1) versus u = aig{fij), 
where fij are numerically obtained largest characteristic multipliers. One can see that with 
increasing delay r, the plot tends to some continuous curve, which determines stability 
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Figure 6: Largest characteristic multipliers for a periodic solution of Duffing oscillator with delay, 
which are rescaled accordingly to the rule 7j = T{\fij\ — 1) (vertical axis) and u! = arg(/Xj) (horizontal 
axis). With the increasing of r, such rescaled spectrum tends the continuous curve 7(cl') representing 
the pseudo-continuous spectrum as described by (f25]l . Parameter values are as in Fig. [3 

of the periodic solution (it is unstable here). Belovi^ we give a proof for the existence of 
pseudo-continuous spectrum. Those readers, who are not interesting in details, may skip 
the following section. 

Proof of the existence of pseudo- continuous spectrum. 

Here we use the theory for linear delay differential equations with periodic coefficients 
, which is analogous to the Floquet theory for ordinary differential equations. This theory 

implies that /i = e^^ is the characteristic multiplier of (l23l) if and only if there is a nonzero 

solution of equation (l23l) of the form 



(26) 



where p{t) = p(t + T) is periodic. Substituting ([26ll into ([23l) . we obtain 

p'it) = {A{t- 1) - Aid) pit) + e-^^B{t- l)p{t - r). 

Since p{t) is T-periodic, we have p{t — r) = p{t — I — nT{l)) = p{t — I) and system (l27l) 
reduces to 

p\t) = {A{t; I) - Aid) p{t) + e-^^B{t; l)p{t - I). (28) 



(27) 



14 



where the large parameter r occurs only as a parameter in e '^'^ . Thus, the corresponding 
monodromy operator of (l28ll 

f/ = f/(A, e"^") 

depends only on A and e"'^'^ smoothly. Since the linear system ( l28ll possesses the periodic 
solution pit) by construction, one of its characteristic multipliers equals to one. This leads 
to the following condition on the monodromy operator 

f/(A, e"^^) - Id] p = 0, 

which must hold for some periodic function p{t). In general case, this is a co-dimension one 
condition, i.e. it leads to some characteristic equation 

F(A,e-^") = (29) 

for determining the characteristic exponents A. This characteristic equation allows proving 
the existence of pseudo-continuous spectrum. Indeed, substituting 

^ = ^ + ^^ (30) 
into (l29ll . we obtain up to the leading order 

F(^|,e-^e-^^-) = 0. (31) 

New unknown real variables now are oo and 7 in Eq. (|3T1) instead of complex A in Eq. 

1 proceed similarly to the case of the pseudo-continuous spectrum for 



( 1291) . Further we wi 
stationary states 
growing phase 



42|. Let us introduce the new artificial parameter instead of rapidly 



F(^^,e-*e-^^) = 0. (32) 



The newly obtained extended equation (l32l) can be generically solved with respect to '^{oj) 
and v^(u;), since the equation is complex, i.e. it contains two real equations for two variables 
7 and (p. The obtained function 7(co') is the resulting asymptotic function, which determines 
the pseudo-continuous spectrum, see Fig. El Coming back from the extended equation f l32] ) 
to the original one ( l3Ti ). we additionally have to take into account the condition 

= p>{u) + 2Txk, A; = 0,±1,±2,... 
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or, equivalently 

- = -(^H + . (33 

IT T 

The last equation (l33l ) determines the discrete values of uj = u^, which correspond to the 
discrete values of the pseudo-continuous spectrum 7(1-1;^). As r increases, the set of ujk covers 
densely the whole domain of u. 

As a result, the pseudo-continuous spectrum approaches the continuous one as r ^ 00. 
Note, that the asymptotically continuous spectrum 7(0;) is determined by a regular system 
of equations fl32|) . which no longer contains the large parameter r. Finally we remark that 
characteristic exponents ( !30ll correspond to characteristic multipliers f l25l) . 



D. Strongly unstable spectrum 

For completeness we show that another type of characteristic multipliers may appear 
which have different asymptotics for large r. These multipliers are not approaching the 
threshold value = 1 as r — > cx), but tending to some unstable value 

Since > 1, the corresponding periodic solution is unstable. In the case, if some periodic 
solution has a strongly unstable multiplier, close-by solutions will diverge on a time interval 
much shorter than the delay r simply because the divergence rate is independent on the 
delay. 

Strongly unstable spectrum may consist of finite number of multipliers (less than n) and 
it occurs if and only if the instantaneous part of ( 123D . i.e. the system of ordinary differential 
equations (ODE) 

i'{t) = A{t-l)m. (34) 

is unstable. The unstable Floquet multipliers of this system will serve as asymptotic values 
Ji for the strongly unstable characteristic multipliers of the system with delay (l23ll . 
An example of strongly unstable spectrum is shown in Fig. [71 

Why does strongly unstable spectrum exist? 

The idea of the proof is the following. Assume that there is a characteristic multiplier ix 
with |/i(T)| > 1, which persists for all r ^ 00 and do not scale with r. Then — > as 
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Figure 7: Example of a strongly unstable spectrum for a periodic solution of Duffing oscillator with 
delay. Single multiplier with large modulus is the strongly unstable one. The others are belonging 
to the pseudo-continuous spectrum, (a) Modulus of largest characteristic multipliers versus branch 
numbers (delay is increasing), (b) Rescaled spectrum, i.e. t(|^| — 1) versus argument of ^ for 
largest characteristic multipliers of a periodic solution on the branch 40. Note the scale on the 
vertical axis. Parameter values: a = 0.5, b = 0.6, d = 0.06. 



T ^ oo and the delayed term in ( 1281 



B{t] l)p{t -l)= ^i'^'^B{t] l)p{t - I) 



(35) 



becomes exponentially small comparing with the term {A{t] I) — \\d)p{t). In fact, choosing 
large enough r, it can be made arbitrary small. Therefore, the equation (l28l) can be formally 
reduced to the ODE 

p\t) = {A{t-l)-\m)p{t). (36) 
The condition for ( l36l) to have a multiplier equal to one reduces to the equation 



det 



t/n - e^^Id 



0, 



(37) 



where Uq is the monodromy matrix of the ODE f l34l ). Since the solutions of (l37ll are multi- 
pliers of ( !34ll . strongly unstable multiplier will approach an unstable multiplier Ji = e^^ of 

m. 



IV. ASYMPTOTIC STABILITY OF BRANCHES 



Let us discuss the main consequences, which follow from the existence of the pseudo- 
continuous spectrum. First of all, let us note, that the sequence of periodic solutions 
x„(t, r) = xoit^l)^ which repeat themselves at time delays / + nT{l), has a well defined 
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1-1 , , , , , , , r- 

1 2 3 4 5 6 7 t 8 



Figure 8: This figure illustrates how branches of periodic solutions have asymptotically well defined 
stability structure. Black solid lines correspond to asymptotically stable parts of the branch for 
large enough delay (—1.22 < ^ < 1.18), dashed are strongly unstable (1.868 < ip < 4.415) gray 
parts are weakly unstable. Parameter values are a = (3 = 2. 

stability limit as n increases. This means, that all solutions from this sequence with large 
enough n will be stable if the corresponding pseudo-continuous spectrum is stable and no 
strongly unstable spectrum is present. Otherwise, the corresponding solutions will be un- 
stable. In other words, in the limit of large delay, the branches of periodic solutions have 
well defined stability structure, i.e. there will be some stable part as well as unstable part. 
The unstable part can be again splitted into strongly unstable (if there are strongly un- 
stable multipliers) or weakly unstable (when only pseudo-continuous spectrum is unstable). 
The corresponding parts can be described by the parameter / on the branches. Figure [8] 
illustrates this using the Stuart Landau model with a = 2 and (3 = 2, see caption to the 
figure. 

Taking into account that almost all parts of the branches are eventually stretching with 
the increasing of delay, an increasing coexistence of stable as well as unstable periodic so- 
lutions is generally expected in systems with large delay. The relative fraction of stable 
solutions depends on a specific system, more exactly, on the asymptotic spectrum distribu- 
tion along the branch. 

V. COMPUTATION OF THE PSEUDO-CONTINUOUS SPECTRUM 

The main equation for finding branches 7(^1^), to which the pseudo-continuous spectrum 
tends to, is given by ( !32i ). As follows from the previous section, the equivalent problem can 
be formulated as follows: 

For any given u, find a value of 7 = 'j{uj) and (p{uj) such that the following linear system 
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with delay 

p'it) = (^Ait; I) - 2^Id) Pit) + e-^-'^B{t; l)p{t - I) (38) 

has a multipher 1. Here I = rmodT and A and B are determined by hnearizing ^ around 
the given periodic solution with period T; see (l24ll . Equivalently, the following extended 
system can be considered 



p\t) = (^D^fixit),xit-l))-^^ld^p{t) 
+e-T-'^D2f{x{t),x{t - l))p{t - I) 



(39) 



with the following additional conditions 



pit) = p{t + T), 
x{t) = x{t + T), 
= 1, 



(40) 



where the first two conditions ensure periodicity and the second one ensures that p{t) is 
nontrivial. This is a typical non-stiff continuation problem, which no longer includes the 
large parameter r. Standard continuation algorithms should be used in order to find the 
functions 'y{uj) and (p^u). 



The implementation of the continuation algorithm will be discussed elsewhere [48||. In- 
stead, we discuss here cases, at which the above problem can be significantly simplified. 

Case 1: rmodT = 0. This situation appears if system ([I]) has a periodic solution at 
r = 0. In this case, the equation (l38ll is reduced to the ODE 



p'it) = A{t; I) - ^-ld + e-T-'^B{t- 1) p{t) 



T 



(41) 



and the equivalent problem for finding functions 'yioo) and (piuj) reduces to the finite- 
dimensional ODE continuation problem [49|, 150||. 

Case 2: system (P) has an additional phase shift symmetry. Examples of such systems 



are Stuart-Landau oscillator (1151 1 or Lang-Kobayashi system [3|, 112| 
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Spectrum 



of some external cavity mode for the Lang-Kabayashi system is shown in Fig. [9l Other 



examples can be found in 



4l|. 
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-1.0 0.0 1.0 -1.0 0.0 irnX 1-0 



Figure 9: Illustration of pseudo-continuous spectrum for external cavity modes of the Lang 
Kobayashi system (adapted from {54 1). 

VI. CONCLUSIONS 



To conclude, we have investigated properties of periodic solutions of systems with delay. 
In particular, we have shown that 

• Periodic solutions of systems with delay are organized in infinitely many branches. 

• The branches of periodic solutions can be obtained as the mapping (JMHI) of a primary 
branch on an appropriate interval of r. From practical points of view, it is sufficient 
to calculate only the primary branch. 

• The branches are eventually becoming wider with the increasing of r, i.e. they occupy 
larger r interval. As a result, the multistability of periodic solutions grows as delay 
increases. 

• This growth of the multistability is linear (fTOl) and the corresponding estimation is 
given in (fTTD and (fT2D. 



One can effectively study asymptotic stability properties of periodic solutions as r 



00. 



As T becomes larger, the spectrum of characteristic multipliers of periodic solutions is 
splitted into two parts: pseudo-continuous spectrum and strongly unstable spectrum. 

The main properties and implications of pseudo-continuous spectrum are explained. 
In particular, pseudo-continuous spectrum controls the destabilization of periodic so- 
lutions. It shows also that the destabilization mechanism of periodic solutions should 
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be similar to that of spatially extended systems . Moreover, it implies that the di- 
mensionality of unstable manifolds of periodic orbits grows linearly as delay increases. 

Strongly unstable spectrum is present when the instantaneous part of the linearization 
around the periodic solution is unstable. In this case, the feedback plays minor role. 



see also 
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In our paper we also outlined the algorithm for practical finding of the pseudo-continuous 
spectrum, although the numerical implementation of the corresponding path-following algo- 



3i- 



rithm (see (l39l ) - (l40l) ) will be discussed elsewhere 
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